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ABSTRACT 



Aims. One of the most challenging and important problem of digital signal processing in Cosmology is the separation of 
foreground contamination from cosmic microwave background (CMB). This problem becomes even more difficult in situations, 
as the CMB polarization observations, where the amount of available "a priori" information is limited. In this case, it is 
necessary to resort to the blind separation methods. One important member of this class is represented by the Independent 
Components Analysis (ICA). In its original formulation, this method has various interesting characteristics, but also some 
limits. One of the most serious is the difficulty to take into account any information available in advance. In particular, ICA is 
not able to exploit the fact that emission of CMB is the same at all the frequencies of observations. Here, we show how to deal 
| with this question. The connection of the proposed methodology with the Internal Linear Composition (ILC) technique is also 
illustrated. 

Methods. A modification of the classic ICA approach is presented and its characteristics are analyzed both analytically and by 
means of numerical experiments. 

Results. The modified version of ICA appears to provide more stable results and of better quality. 
Key words. Methods: data analysis - Methods: statistical - Cosmology: cosmic microwave background 



1. INTRODUCTION 



The experimental progresses in the detection of cosmolog- 
ical emissions require a parallel development of data anal- 
ysis techniques in order to extract the maximum physical 
information from data. In particular, different emission 
mechanisms are characterized by markedly distinct un- 
derlying physical processes. Data analysis often requires 
the component separation in order to study the individ- 
ual characteristics. To achieve such a goal, a link between 
the branch of signal processing science which character- 
izes and separates different signals and astrophysics is yet 
well established, and in many cases, modern signal pro- 
cessing techniques have been imported and applied in an 
astrophysical context. This is the case of the Independent 
Component Analysis (ICA), used for the separation of the 
Cosmic Microwave Background (CMB) fro m diffuse fore- 



ration of all the components contributing to an observed 
signal. 

More formally, let the available data in the form of N 
mean-subtracted maps {y}fLi, corresponding to different 
observing channels, and containing M pixels each. If Y i = 
VEC T [^i] 0, these maps can be arranged in a N x M 
matrix 



Y 2 



\Y N J 



(1) 



A common assumption in CMB observations is that each 
yi is given by the linear mixture (or in astrophysical 
term "frequency channel") of N c components {«S J }^ 1 
due to different physical processes (e.g., free-free, dust re- 



grounds originated by our own Galaxy (see IStivoli et al. 



20061 and references therein). This techniques offers sev- 
eral advantages. In particular, under the only assumption 
of mutual statistical independence, it permits the sepa- 



1 We recall that the operator VEC[i/] transforms a matrix 
H into a vector by stacking its columns one underneath the 
other. Hence, VEC T [H] provides a row array. 
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radiation, bremsstrahlung, 



In formula, 



with Sj = VEC [Sj] and {<%•} constant coefficients. 
With this model, it is hypothesized that for the jth phys- 
ical process a template Sj exists that is independent of 
the observing channel " i ". Although rather strong, ac- 
tually it is not unrealistic to assume that this condition is 
satisfied when small patches of the sky are considered. In 
matrix notation, Eq. ([2]) can be written in the form 



where 



Y = AS, 



( Si \ 

s 2 



(4) 



and 



\SnJ 
( an a 12 ... a 1Nc \ 

&21 CL22 • • • &2N C 
\ &N1 CLN2 • • • &NN C ) 



(5) 



denotes the so called mixing matrix. 

Here the problem is that only the mixtures Y are avail- 
able, whereas neither A nor S are known. Hence, the issue 
raises whether from Y it is possible to obtain the compo- 
nents {Si}. Surprisingly, a positive answer is possible. At 
this point, however, it is necessary to stress that the prob- 
lem presents a basic ambiguity. In particular, at best each 
Si can be determined unless a multiplicative constant. In 
fact, if Si is multiplied by a scalar and the corresponding 
i-th column of A is divided by same quantity, then an 
identical model is obtained. For this reason, it is custom- 
ary to assume that the variance of Si is equal to one, i.e., 
SiSj/M = 1. 

For simplicity, in the following it is assumed that N c = 
TV, i.e. that the number of observed mixtures is equal to 
the number of components. In this way, A is a square 
matrix. 



2. A CLASSIC APPROACH: ICA 

One of the most celebrated technique for the blind sep- 
aration of signals in mixtures is the so called indepen- 
dent component analysis (ICA). The basic idea behind 
ICA is rather simple (and obvious): to obtain the sepa- 
ration of the components it is sufficient to have an es- 
timate of A' 1 . In fact, S = A~ X Y . Now, if the CMB 
component and the Galactic ones are mutually uncorre- 
cted, i.e. if the corresponding covariance matrix is given 
by C s = SS T /M = I, then 



This system of equations defines A but an orthogonal 
matrix. In fact, if A = ZV, with V orthogonal, then 



YY 1 /M = AA 



ZVV T Z T 



ZZ . The problem is 



ft' that, given the symmetry of YY T , system © contains 
only N(N + l)/2 independent equations, but the esti- 
mates of N 2 quantities should be necessary. In ICA, the 
N (N — l)/2 remaining equations are obtained by enforc- 
ing the constraint that the components {Si} are not only 
mutually uncorrelated but also mutually independent. In 
other words, the separation problem is converted into the 
formal 

S = argminF(S'), (7) 
s 



(3) subject to S = A~ X Y and YY T /M = AA T , (8) 



with F(S) a function that measures the independence be- 
tween the components {Si}. The definition of a reliable 
measure F(.) is not a t rivial task. In literature , various 
choices are available (see iHvvarinnen et af. 2001, and ref- 
erence therein). In practical algorithms, the optimization 
problem is not implemented explicitly in the form ([7|)-(j8j). 
Typically, a first estimate S* is obtained through a prin- 
cipal component analysis (PCA) step followed by a spher- 

ing operation (i.e. forcing S*S^ /M = 1). In this way, 
a set of uncorrelated components become available with 
Cg = /, as well the corresponding mixing matrix A*. 
Later, these estimated quantities are iteratively refined to 
maximize F(S) and to get the final S. Again, in litera- 
ture, various techniques are available. Among these, one 
of the most famous and used algorithm i s FAS TIC A based 



on & fixed-point opt i mization approach (IHvvarinnen et al. 



200l[ lMainoll200l iBaccigalupi et al. 112004 ). An alterna- 
tive technique is JAD E that makes u se of the joint diag- 
onalization algorithm (jCardoso II1999I ). 

3. A SUBSPACE APPROACH 

The main benefit of ICA is the fact that, apart from the 
independence of the components {Si}, it does not make 
use of further assumptions. If from one side, this makes 
the method easy to use, on the other one it does not per- 
mit the exploitation of any information that be available 
in advance. In particular, in the case of the CMB observa- 
tions, it is expected that, at least on small patches of the 
sky, the mixing matrix A can be written in form: 



A = 



(I 
1 



d\2 

0,22 



aiN \ 

d2N 



(9) 



\ 1 &N2 • • • &NN J 



This means that the component S\, here assumed to cor- 
respond to the CMB emission, gives the same contribution 
in all the observed mixtures. Here, the question is how to 



YY T /M = AA 7 



2 We recall that the function "arg min F(x)" provides the 
values of x of for which the function F(x) has the smallest 



(6) value. 
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implement this piece of information. A possible solution 
can be obtained if model (|3j) is written in the form: 





/1\ 




( a 12 \ 




1 






Y = 




Si + 










\ a N2 J 



( a 1N \ 

\ Q>NN / 



(10) 



This equation enlightens the fact that the columns of A 
span a signal space < A >= span{l, 02, . . . , ajv} where 
1 = (1,1,..., 1) T and a,i denotes the i-th column of A. 
In other words, signal Y lives in a iV-dimensional space. 
Moreover, from the same equation it is evident that if Y is 
projected onto a (TV — l)-dimensional subspace orthogonal 
to vector 1, then the contribution of £1 is removed from 
Y itself. This job can be done by means of the projection 
matrix 

P_L = I- 1(1 T 1)- 1 1 T . (11) 

The application of ICA to 



Y± = P±Y, 



(12) 



permits to obtain the estimates £2 , £3 , • • • ,Sn of the 
N — 1 components S 2 ,Ss,--- ,Sn- Hence, since it is 

1 A — Sij (Sij denotes the Kronecker func- 

tion) and by assumption {SiS i }fL 2 = 0, the columns of 
A can be estimated through Eq. ([TO]) by means of 



{S i S j /M}i J=2 



ai = YSt /M, 



2,3,. ..,7V. 



(13) 



At this point, always using Eq. ([TO]) . £1 can be derived 
from 

Si = 1 T (Y -A^S-^/N, 



where 
and 




Since SiS i = 0, i ^ 1, it is not difficult to see that 
Si = Si- 



(14) 
(15) 

(16) 
(17) 



It is worth noticing that, contrary to the other compo- 

nents, the variance of £1 is not one: S\S 1 /M ^ 1. This 
is a consequence of the fact that a\ = 1 is fixed. 

Here, it is necessary to stress that, if one is interested 
only in the component Si, then the situation is simpler 
since it is not necessary that the other components are 
independent of even uncorrelated. The point is that the 
separation of the components {Si}fL 2 has no effect on 
the computation of Si. In other words, it does no mat- 
ter whether these components are correctly disentangled 
or not. In fact, the same £1 as in Eq. ([Tl]) is obtained if, 
instead of the independent {Si}^L 2 ^ m Eq. ([13]) . the or- 
thogonal {S^ }^ 2 are used that are computed through the 



application of the PC A to Y±. This is because Eq. ([TO]) 
requires the orthogonality of the components, not their 
independence. Moreover, Eq. ([3]) can be written in the 
equivalent form 

Y = AS, (18) 



where 



A = AH, 

... 0\ 



H 



/I 





(19) 
(20) 

(21) 



and $ is any arbitrary (N — 1) x (TV — 1) non-singular 
matrix. Now, since Si = #1 and the first colum of A is 
still 1, the meaning of this equation is that there is an 
infinite number of sets {^1^2 that, when projected onto 
the N — 1 subspace orthogonal to 1, produce the same 
Y±. Hence, for the separation of Si it is not necessary 
the use of the true {S^}^ but only one of such sets. 

Of course, these are theoretical results. In practical sit- 
uations, it is quite improbable that for finite signals the 
condition {SiSj}i = 0, i 7^ 1 be strictly satisfied. In fact, 
because of the statistical fluctuations, in general the com- 
ponents {Si}^ present a certain degree of mutual corre- 
lation even in the case they are the realization of indepen- 
dent, stationary, stochastic random processes. Therefore, 
enforcing the condition SiSi = 0, i ^ 1, makes inaccu- 
rate the estimate of the coefficients {a^} as provided by 
Eq. HH 



4. SOME NUMERICAL EXPERIMENTS 

Because of the arguments presented above, it is to be ex- 
pected that, with respect to the classic ICA, the use of 
a subspace method provides more accurate and stable re- 
sults. Given the non-linear nature of the algorithms, the 
verification of this expectation has to be made through 
numerical experiments. 

Here, non- astronomical subjects have been deliber- 
ately chosen. In this way, a direct visualization of the sepa- 
ration is possible and hence an easier and safer assessment 
of its goodness. Moreover, the use of deterministic subjects 
make easier the modeling of various experimental condi- 
tions (e.g., the sample correlation between different images 
can be obtained using images with almost-constant lumi- 
nosity areas in correspondence to the same coordinates). 

Figures [l][2] show the results provided by ICA and its 
version based on the subspace approach (modified ICA) 
when three mixtures are available each containing the 
contribution of an equal number of (almost) uncorrelated 
components with a mixing matrix 



1.0 0.2 0.4 
1.0 0.5 0.3 
1.0 0.3 0.6 



(22) 
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No noise has been added. The examination of these figures 
seems to indicate that both techniques are able to produce 
an excellent separation. Actually, the results provided by 
the classic ICA are unstable and often unsatisfactory. The 
point is that this method is non-linear and therefore the 
algorithms have to be initialized. As a consequence, differ- 
ent results can be obtained according to the chosen initial- 
ization. This is quantified in Fig. 02 where the normalized 
norm of the residuals, \\Si — Si\\/\\Si\\, corresponding to 
200 different initializations are shown. The residuals are 
computed from the difference between the real solution 
and the estimated one. From this figure, it is clear that 
the classical ICA provides more more than one solution, 
while the modified ICA method produces a stable Si. This 
is a consequence of the fact that in the subspace approach 
such component is computed via linear operations only. 

The question of the stability of the solution may be 
considered of secondary importance. Actually, this is not 
true. The ICA method does not permit to check if a spe- 
cific solution is satisfactory or not. Any separation pro- 
vides components that, when summed up, will perfectly 
reproduce the original mixtures. In a real experimental 
situation, this means the unavailability of a reliable se- 
lection criterion. For example, a simple selection criterion 
could be based on the frequency with which a solution 
is obtained for different initializations of the algorithm. 
However, there is no guarantee that the most frequent so- 
lution represents the best one. This is the case for the 
experiment in Figs. 0H6] where, contrary to the previous 
one, the components present a certain degree of correla- 
tion with a (normalized) cross-product matrix 



C s 



1.00 0.20 0.06 
0.20 1.00 0.06 
0.06 0.06 1.00 



(23) 



As expected, the separation is by far less satisfactory than 
in the previous experiment. However, also in this case the 
results concerning the subspace method appear more sta- 
ble. Here, the point of interest is that from the exami- 
nation of the top-left panel in Fig. [6] it is possible to see 
that the classic ICA provides the estimate Si closest to 
the true solution. However, apart from the fact that the 
#2 corresponding to such estimate is systematically the 
worst among those obtained, the frequency with which the 
"bestf Si is found is very small. In a practical situation, 
such a solution should have been discarded. 

It has not to surprise that, when the conditions of ap- 
plicability are violated, the classic ICA method is able to 
"see" a good estimate of Si that is "unreachable" with the 
subspace approach. This is because, as stated above, the 
classic ICA works with a larger number of unknowns and 
therefore is more flexible and can span a wider "solution 
space" . However, there is no guarantee that such a flexibil- 
ity is effectively fruitful. The situation is similar (although 
not identical) to a polynomial fit: the use of high degree 
functions permits a greater flexibility but at the cost of a 
remarkable instability of the results that often make quite 
hard, if not impossible, the choice of a good solution. 



5. RELATIONSHIP WITH THE ILC METHOD 

In CMB literature, another method has been often used 
for the separation of Cosmic signals from the Galactic fore- 
grounds. This i s the so called internal linear combination 
method (ILC) teejmetlel_ajj l20oj; lEriksen et aTll20Q4l : 
Iffinshaw et al. ll2QQ7h . The aim of this method is not the 
separation of all the signals that contribute to the observed 
mixtures, but only the extraction of the specific compo- 
nent Si. With ILC a solution is searched in the form 



Si = w Y, 



(24) 



with the column vector w providing a set of appropriate 
weights. Since, the basic assumption is that Si is the same 
in all the mixtures, i.e., 



Y = l^i+ AT; 



(25) 



with Af a zero-mean noise that provides the contribution 
of all the components other than that of interest, it is 
imposed that 

l T w = 1. (26) 



In this way, 



Si = Si + w T Af, 



(27) 



i.e. the weights do not alter the Si component. For the 
same assumption, among all the possible solutions pro- 
vided by Eq. ([24|) with the condition ([26|) , that of interest 

has the property that a 2 = SiS 1 is a minimum. In fact, 
assuming the noise Af uncorrelated with S, it is 



1 



a 2 = —[SiS{+w T ATAr T w]. 



M 



(28) 



Hence, the minimization of a 2 with respect to w implies 
the strongest filtering of the component Af. It can be 
shown th at the weights which minimize this quantity are 
given by (|Eriksen et al. 1120041 ) 



w 



(29) 



where Cy = YY T /M. Hence, the ILC estimator takes 
the form 



Si = 



\ t Cy X Y 



(30) 



Although this estimator appears different from that given 
by Eq. ([T4]) . actually they provide identical results. In fact, 
under model (HOI), it is 



CV = ACsA 1 



If this equation is inserted in Eq. (|30|) . one obtains 

Si = al T A^C^S, 
with the scalar a given by 

a=[l T A- T Cs 1 A- 1 l}-\ 



(31) 



(32) 



(33) 
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and A T = (A 1 ) T . Now, since it is trivially verified that 
1 T = ej A T , 

where 

ei = (1,0,...,0) T ) 
it is 

1 T A~ T = ej. (36) 

Hence, a = (Cg 1 )n = ([SiSj]) -1 = cr^. As a conse- 
quence, if Si is uncorrelated with {Si}^L 2 i.e. if 



(34) 
(35) 



cm 





. 


• ) 


the form 





0~22 


0"23 • 




' (37) h A • 

where A is 





0~N2 


0~N3 • 


. (JNN ) 


Ai > A 2 > 



one side, the case N < N c does not offer many possibil- 
ities to obtain meaningful results, on the other one the 
case N > N c can be successfully addressed. In fact, N c 
can be determined by the number of non-zero eigenvalues 
of matrix Cy and the corresponding eigenvectors can be 
used to construct a basis of the signal-space with the cor- 
rect dimensionality. After that, it is sufficient to project 
Y onto this space, obtaining a "new" set of N* < N mix- 
tures Yr, and then to work with this. The rest of the 
procedure in Sect. [3] and Sec. [5] remains the same. More 
in particular, if the N x N matrix Cy is decomposed in 



C s 



from Eq. ([32]) and the fact that C s l has the form 



VAV T 



(42) 



V o 







(38) 
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with fl the bottom-right block of the matrix in the rhs of 
Eq. (j37J), one obtains that 



S 1 = efS = S 1 , 



(39) 



i.e., the same result as Eq. (QZj). More in general, the es- 
timators ([HI) and ([32]) provide identical results also when 
the components Si is not uncorrelated with the other ones 
and/or instrumental noise is added to the observed mix- 
tures. The reason is that, as stated earlier, ILC provides 
an estimate Si with the property that the quantity a 2 
in Eq. ([28]) is a minimum. Although not evident from 
the treatment in Sec. [3l the same holds for the estima- 
tor (QZj). In fact, it is not difficult to realize that, after the 
determination of the components {Si}fL 2 , the coefficients 
{a^}^ 2 , as given by Eq. ([T3]h are the solution of 



da 2 
dA_ 



dAi 



= 



with 



(40) 



(41) 



SiS, =1 T (Y- A_i£_i) T l 
that is derived from the sample version of Eq. ([TQ|) . 

6. FINAL REMARKS 



In the previous section it has been assumed that the num- 
ber N of the observed mixtures (or frequency channels, 
images at different observing frequency) equals the num- 
ber N c of the components. In practical application this co- 
incidence is improbable. Of course, in order the subspace 
approach can work satisfactorily, it is necessary to know 
the correct dimension of the signal-space. Therefore, the 
question raises on what happens when A^ c =^ N. If from 



> ^n c > ^n c +i = • • • = Ayv = 0, whereas 
V is an orthogonal matrix whose columns contain the cor- 
responding eigenvectors, then 



R 



: 1/2 v t y, 



(43) 



where A 1 ^ 2 is a diagonal matrix whose non-zero values 

1/2 

are given by the inverse of the non-zero entries of A 1 . 

More complex is the situation when Y is contaminated 
by measurements errors N. If N is zero- mean and addi- 
tive, model ([3|) converts into 



Y = AS + N. 



(44) 



Here the problem is that the effective number of com- 
ponents becomes N c + N and therefore the separation 
problem is always underdetermined. Therefore, although 
E[N] = 0, it h appens that E\S}] ^ S \. In other words, a 
bias is present ( Vio fc Andreanill2QQ8l ). This point is more 
evident if the ILC estimator ([30]) is considered. There, the 
bias is due to the fact that matrix C y is no longer given 
by Eq. but by 



Cy = AC sA T 



' N • 



(45) 



This equation suggests that a simple way to remove the 
bias is to use Cy — Cat instead of Cy • Something similar 
holds also for S\ as provided by Eq. ([T4|) . In fact, after 
some algebra it is possible to show that an equivalent form 
is 

Si = 1 T [I - CyP1(P±C y p1?P±]Y/N, (46) 

where symbol " t " denotes pseudo-inverse. Hence, the 
same arguments apply as above. Concerning the influ- 
ence of the noise of the components {Si}, i / 1, the 
situation is much more difficult since similar to that en- 
countered in the classic ICA approach (for details, see 
Hvvarinnen et al. Il200ll ). 



7. SUMMARY 

In this paper we have considered the problem of a modi- 
fication of the ICA separation technique that permits to 
exploit the u a priori' information the, contrary to the 
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Galactic components, the contribution of CMB to the mi- 
crowave maps is independent of the observing frequency. 
A subspace approach has been proposed that is more sta- 
ble and provide more accurate results than the classic ICA 
technique. A relationship between this approach and the 
Internal Linear Composition method has been also shown. 

References 

Baccigalupi, C, et al. 2004, MNRAS, 354, 55 
Bennett, C.L., et al. 2003, ApJS, 148, 97 
Cardoso, J.F. 1999, Neural Computation, 11, 157 
Eriksen, H.K., Banday, A.J., Gorski, K.M., & Lilje, P.B. 

2004, ApJ, 612, 633 
Hinshaw, G., et al. 2007, ApJS, 170, 288 
Hyvaxinnen, A., Karhunen, J., & Oja, E. 2001, 

Independent Component Analysis (New York: John 

Wiley & Sons) 
Maino, D., et al. 2002, MNRAS, 334, 53 
Stivoli, F., Baccigalupi, C, Maino, D., & Stompor, R. 

2006, MNRAS372, 615 
Vio, R., & Andreani, P. 2008, A&A, submitted 



R. Vio, & P. Andreani: ICA method 7 

Original N.1 Original N.2 Original N.3 




50 100 150 200 250 50 100 150 200 250 50 100 150 200 250 



Observed N.1 Observed N.2 Observed N.3 




50 100 150 200 250 50 100 150 200 250 50 100 150 200 250 



Fig. 1. Top panels - Original images S used in the experiment dealing with component separation as described in 
Sect. HJ Bottom panels - Observed mixtures Y = AS, with A given by Eq. ([22]) . In this experiment the components 
are almost uncorrelated. 
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Reconstruction (ICA) Reconstruction (ICA) Reconstruction (ICA) 




50 100 150 200 250 50 100 150 200 250 50 100 150 200 250 



Reconstruction (Modified ICA) Reconstruction (Modified ICA) Reconstruction (Modified ICA) 




'V- ^ A , i J ^DU I" ■ ^{J im:.x.x.-^>.o>x^o^> : o>....,>:.:,.>x^..o.c,v.: v :.xoxo> : ^ 

50 100 150 200 250 50 100 150 200 250 50 100 150 200 250 



Fig. 2. Top panels - Typical separation obtained with the classic ICA algorithm when applied to the observed 
mixtures shown in Fig. [TJ Bottom panels - Typical separation obtained with the subspace based ICA when applied 
to the observed mixtures shown in Fig. [U 
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Classic ICA 



Modified ICA 
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300 




150 



100 



150 



100 



7.45 



x 10" 




x 10" 




0.5 1 

Normalized Norm of the Residua^ Q -3 



Fig. 3. Histogram of the normalised residuals \\Si — The residuals are defined as the difference between 

the real and the estimated solution. The corresponding solutions of the component separations are obtained with 200 
different inizializations of the classic ICA and subspace based ICA methods when applied to the observed mixtures in 
Fig. fTJ top panels correspond to component i = 1, central panels to i = 2, bottom panels to i = 3. As shown in the 
top panel the solution in the modified ICA approach for component i = 1 is stable. 
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Original N.1 



Original N.2 



Original N.3 






50 100 150 200 250 



50 100 150 200 250 



50 100 150 200 250 



Observed N.1 



Observed N.2 



Observed N.3 




50 100 150 200 250 




50 100 150 200 250 




50 100 150 200 250 



Fig. 4. Top panels - Original images S used in the experiment dealiong with component separation as described in 
Sect. IU Bottom panels - Observed mixtures Y = AS. In this experiment the components are partially correlated 
with the correlation matrix given by Eq. ([23]) . 
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Reconstruction (ICA) Reconstruction (ICA) Reconstruction (ICA) 




50 100 150 200 250 50 100 150 200 250 50 100 150 200 250 



Reconstruction (Modified ICA) Reconstruction (Modified ICA) Reconstruction (Modified ICA) 




50 100 150 200 250 50 100 150 200 250 50 100 150 200 250 



Fig. 5. Top panels - Typical separation obtained with the classic ICA algorithm when applied to the observed 
mixtures in Fig. [4j Bottom panels - Typical separation obtained with the subspace based ICA when applied to the 
observed mixtures in Fig. [H 
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Fig. 6. Same as Figure [2] for the simulations shown in Figure [H and [51 



